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Abstract. The parareal in time algorithm allows to perform parallel simulations of time depen- 
dent problems. This algorithm has been implemented on many types of time dependent problems 
with some success. Recent contributions have allowed to extend the domain of application of the 
parareal in time algorithm so as to handle long time simulations of Hamiltonian systems. This 
improvement has managed to avoid the fatal large lack of accuracy of the plain parareal in time algo- 
rithm consequence of the fact that the plain parareal in time algorithm does not conserve invariants. 
A somehow similar difficulty occurs for problems where the solution lacks regularity, either initially 
or in the evolution, like for the solution to hyperbolic system of conservation laws. In this paper 
we identify the problem of lack of stability of the parareal in time algorithm and propose a simple 
way to cure it. The new method is used to solve a linear wave equation and a non linear Burger's 
equation, the results illustrate the stability of this variant of the parareal in time algorithm. 

Key words, parareal in time algorithm, parallelisation, time discretization, evolution equations, 
hyperbolic system, wave equation. 

AMS subject classifications. 65L05, 65P10, 65Y05 

1. Introduction. Parallel in time algorithms represent a new competitive way 
of using the ever increasing number of cores in today's supercomputers platforms in 
cases where classical domain decomposition methods either are of no relevance (for 
systems of differential equations) or in cases where the scalability of the sole domain 
decomposition (for evolution PDE's) gets to saturation. This kind of algorithm in 
the 4 th direction has received too little attention in the literature in the past due 
to the inherent sequential nature of the time direction: it seems indeed difficult to 
simulate what will arrive in the far future, e.g. next week, without knowing in details 
what will arrive in the close future, e.g. tomorrow and each day of the current 
week. Due to the paramount importance of finding new ways for parallelisation, 
going together with the ever increasing number of processors present in the new 
generations of architectures, algorithms such as parareal in time, parallel deferred 
corrections or parallel exponential integrators have become an active research topic 
showing promising applicability for large-scale computations and perhaps one way 
to go to exascale computations. We refer e.g. to the book of K. Burrage [I] for a 
synthetic approach on the subject (see also [5]). In this book, the various techniques 
of parallel in time algorithms are classified into three categories: (i) parallelism across 
the system, (ii) parallelism across the method and (iii) parallelism across the time. 

The parareal in time algorithm pioneered in [M] extended later in [H El [18] under 
a form much better tuned to the treatment of nonlinear problems enters in the third 
category above (see [3U] for a first algorithm in this category). It combines two 
propagators, one fine and one coarse, in a predictor-corrector manner that allows 
to use the fine propagator in parallel on the various processors, while the coarse 
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propagator is solved in a sequential manner. The algorithm provides good, provable 
convergence, for diffusive PDEs such as parabolic type problem. However, several 
studies have shown some instabilities for hyperbolic equations [51 [IT] [Jj3 [5] and the 
numerical analysis performed on simple examples or partial differential equations in 
[3j [2j [21] allows us to understand this behavior. 

A challenging field of application for parallel in time propagators is for long time 
simulation of Hamiltonian systems. Here also the plain method lacks the geomet- 
ric properties that are essential to guarantee the quality of the approximation. Let 
us remind that these geometric properties lead to the convergence of the algorithm 
towards the solution of the original systems on long time simulations. These proper- 
ties — either the symplecticity or the symmetry of the scheme — are not preserved 
for the parareal scheme even if the basics coarse and fine propagators are symplectic 
or symmetric. A cure to this lack of geometric properties has been proposed in [TJ, 
based on two ingredients: the first one is the symmctrization of the parareal in time 
algorithm that leads to a new multistep schemes. This ingredient alone is not suffi- 
cient though to provide a correct algorithm for Hamiltonian systems. Indeed, some 
resonances are artificially introduced by the parareal strategy itself as is explained 
in [TJ; they prevent the symmetric variant of behaving well. The second ingredient 
that has been introduced in [7J consists simply, after each correction of the predictor 
corrector scheme, to project the solution coming out the parareal in time algorithm, 
on the manifold expressing the preservation of some basic invariant quantities of the 
Hamiltonian system. This is quite simple to implement, at least in its basic formula- 
tion (there exists a symmetric version based on a more involved implicit resolution of 
projection on the manifold). 

Following this second idea advocated to stabilize the parareal in time algorithm 
for Hamiltonian system, we propose in this paper an extension of it allowing to cure 
the lack of stability that occurs when solving hyperbolic type problem. Another 
strategy, based on providing a much more involved correction procedure re-utilizing 
previously computed information has been proposed in [9] [10] and already improved 
the performance of the parareal in time algorithm for second-order hyperbolic systems. 
The method proposed here seems much more simple to implement and more general 
in scope and applications. 

The remainder of this paper is organized as follows. In Section 2, we remind the 
basics on the plain parareal in time algorithm including the current state of the art 
on the stability and convergence analysis of the algorithm; we use this to introduce 
the proposed stable variant of the parareal scheme. In Section 3 we address the case 
of the wave equation, a simple second-order hyperbolic equation. The correction we 
propose is discussed in full details and the numerical simulation illustrates the good 
properties of the variant. In Section 4 we address the case of the Burgers' equation, 
with small viscosity, here again the variant performs well. In Section 5 we present 
some elements of numerical analysis proving the stability of the proposed variant of 
the parareal scheme and finally some conclusions are drawn in Section 6. 

2. Some elements on the parareal in time algorithm. 

2.1. Basics on the parareal in time algorithm. Let us consider a possibly 
non-linear time dependent partial differential equation of the form 



-£+A(t;y) = 0, y(0) = y 



(2.1) 



Stable parareal in time algorithm for hyperbolic system 



3 



where y(t) takes its values in a Banach space B and A is a partial differential operator 
defined over B. We assume here that this problem is well posed, at least for some 
time: the solution y(t) exists for t G [0,T). More precisely, we assume that there 
exists a propagator £ such that 

y(t) = elivo), 

named exact propagator. We note that, with obvious notations, we have 

y(t) = £l(y(s)), Vs,t, < s < t < T. (2.2) 

Let us assume now that we want to simulate this PDE on the interval [0,T), 
we thus introduce a discrete propagator J- that is assumed to be a fine and precise 
approximation of the exact solution based on a discretization in time only. The 
propagator T can be based on an approximation using a time discretization and, e.g. 
the use of an Euler-type or a more involved scheme based on a small enough time 
step St. 

For instance, let us denote as T = 0, T N = T and T n = nAT (with AT = ^) 
special times at which we are interested to consider snapshots of the solution y(T n ), 
then, with notations similar as above (|2.2[) . we have the approximation 

Y n =J r To(yo), n=l,...,N 

we also have 

Yn+l = J~^ + 1 {Y n ). 

In what follows we present the parareal in time algorithm that allows to define a 
sequence {Y k } k ,Y k = (Yf , Y%, Y$), such that, for any n, Y k converges to Y„ in 
the space B as k goes to infinity. The definition of this sequence requires the avail- 
ability of another approximated propagator, denoted as Q. This second propagator 
is assumed to be cheaper than T at the price of some inaccuracy. One can think im- 
mediately about having Q based on a similar Euler scheme as T with the larger time 
step dT, St « dT < AT. But other possibilities are also offered: in addition to the 
solution scheme in time, a discretization in space has to be performed, for instance, 
the spacial discretization can be based on a finite element method, the coarse solver Q 
can then be based both on a coarser time step and a coarser spacial discretization, we 
refer to [TT] and [T7] where such a strategy is followed. Another possibility consists 
in providing J- with all the physics of the phenomenon while Q will be based on a 
simplified physics (see eg [16] for an example). 

Assume T and Q are given, we can proceed to the definition of the sequence. Let 
us initiate Y° by 

Y^=g^(y ), n = l,...,N (2.3) 
Assume now Y fe is known, the iterative process proceeds by setting Yq = y(0), and, 

= Gr: +1 (Yn +1 ) + FtT^) ~ Gl: +l {YZ)> n=l, ...,N. (2.4) 

We first comment on the parallelization of the above algorithm. The initialization 
step (|2.3|) is sequential but involves the cheap coarse propagator. Then assuming 
that every (Y k ) n is known at step fc, we first solve over each ]T n ,T n+ i[ the totally 
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dt 

r i i 

dT 



independent propagation problems that provide J-" T ™ +1 (Y*); we can either propagate 

again (and now in parallel) the coarse solver to get Q^ n+1 (Y*), or just remember that 
these quantities were computed during the previous step and stored, any way this is 
again parallel work. The corrector (Y*) — G^ n+1 (Y*) can then be built, for any 

n and the step can proceed. We finally solve the predictor based on the coarse solver 
Gt™ +1 (Yn +1 ) — sequentially — and immediately correct it by adding the previously 
computed corrections until convergence. Note that in this process, the fine solutions 
are only used on small propagation intervals and can be performed in parallel. 

2.2. Numerical analysis of the parareal in time algorithm. We review 
in this subsection some elements on the stability and convergence of the parareal in 
time algorithm that have been first presented in [2]. Let us assume that there exists 
a sequence of Banach spaces Bj such that -Bj+i C Bj. It is also assumed that the 
problem I2.1l is stable in these spaces, in the sense that there exists a constant C > 0, 
such that, for all j,l < j < J 

\\%(Vo)\\ Bi < C\\yo\\ Bj (2.5) 

The algorithm converges, under the following natural hypothesis: denoting by SJ- 
(resp. 5Q) the difference S T l s — £* — J 7 * (resp. SQ* S — £ l s — (/*), we assume that there 
exists a constant C > 0, such that, for all j, 1 < j < J 

Vi > 0,Vr > 0, \\gt +r (x) - Qt+ r (y)\\ Bj < (1 + Ct){\\x - y\\ Bj ). (2.6) 

In addition let us assume that there exists a constant C > 0, such that, for all 

i, i < j < J 

Vi>0,Vr>0, \\5F t t +T (x)-5T t t +T { V )\\ B] <CTi 1 \\x-y\\ B]+ll (2.7) 
\\SQ t t +T (x)-Sg t + T (y)\\ Bj < Cre\\x-y\\ Bj+1 , (2.8) 

where r\ and e are small quantities, typically, for the Euler scheme these are r\ = 6t 
and e — AT and C is a generic constant, independent of the quantities on the right 
hand side of the previous inequalities and on j < J, but may depend on J. 

Remark 1. For Euler schemes based on a time step St for J- and AT Q , hypoth- 
esis {2.7}) and (2.8]) hold with rj = St and e = AT respectively. 

The convergence is stated in the following (see theorem 1 in [2]) 
Theorem 2.1. Assume that the discrete propagators J- and Q satisfy \2.6}) . {2.7}) 
and i2.8}) . then the error between the exact solution and the solution provided by the 
parareal scheme {2.4\ l satisfies for any k > 0: 

Vn<N, ||y„ fc -y(r")|| Bo <C( £ fe+1 +r,)(l + ||y || Sfe ). (2.9) 
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Proof: The proof is done by induction over j > for the statement 

Vi>0, \\Y4-y(T n )\\ Bi <C(^+ 1 +T,){l + \\yo\\ Bi+j+1 ), (2.10) 

which is obvious for j = since it corresponds to the analysis of the plain coarse 
operator. Let us assume that it is true for j and let us prove it for j + 1. From the 
definition of Y^\ , we write 

= Ql: +1 {Yt l )-QTl +1 {y{T n )) 
+[s^ - g^m) - - gl: +1 ]{y{T n )) 

-[£^ + - j£ +1 ](y(T")) 

or again with the notations introduced above 

Yiti - y(T n+i ) = - gi: +i {y{T n )) 

+sg^ 1 (Y^)-sg^ 1 ( y (T n )) 

-S^: +1 (Y r i)+6^: +1 (y(T n )) 
-[4: +1 -^; +1 ](y(T")) 

hence by triangular inequality 

\\Yj+l - y(T" +1 )|| Bi < \\g% +1 (Y-> +1 ) - gr: +1 (y(T n ))U 
+ \\SgT: +1 (Y^)-5g^(y(T n ))\\ Bt 
+ \\6^(Yi) + 6^(y(T n ))\\ Bi 

+\\[£tT - ^: +1 ]{y{T n ))\\ B , 

that from hypothesis (|2.6[) . (|2.7I) and (12.81) can be upper bounded as follows 

\\Y^l-y(T n+1 )\\ Bi < (1 + CAT)\\Y,{ +1 -2/Cn|U 
+CAT £ ||^-y(T")|| Bi+1 

+CAT»,||^ -J/(T»)|| Bi+1 +CAT7,||2,(T n )|| Bi+1 

the stability hypothesis (|2.5p together with the induction hypothesis (|2.10[) allow to 
get 

\\Yitl ~ y(T n+1 )\\ Bl < (1 + CATWt 1 - y(T n )\\ Bi 

+CAT(e + T))(eJ+ l + v)ho\\B l+]+2 + CAT V \\y \\ Bi+1 
< (l + CATWZ+'-yiTnWB, 
+CAT(ei+ 2 + V )\\y \\B 1+J+2 

from which we get (|2.10l) for j + 1 from the discrete Gronwall lemma, and the proof 
is complete. 
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It is interesting to notice two things at this level. First the constant C in (|2.9p 
depends on k and T and this convergence result can be polluted on a long time range 
by a constant C that becomes too large. Second, the convergence results implies, 
of course, stability. Hence, provided that the initial condition is regular enough and 
remains regular the parareal in time algorithm is stable. 

The previous remark explains actually the problem of the use of the parareal 
scheme when non dissipative PDE or PDE with small dissipation are simulated. This 
is the case for small viscosity parabolic problem as is explained in the next subsection. 

This is also the case for hyperbolic systems. We can illustrate this statement 
on two equations: the periodic wave equation and the periodic Burger's equation for 
which, starting with regular enough initial condition solves or corrects the instability 
of the parareal in time algorithm in this situation. This allows to rectify somehow 
the statements that were done before in e.g. [TU] where the stability problem occurs 
for second-order hyperbolic problems and not for parabolic and first-order hyperbolic 
problems. The problem is not really on the type of the equations (and certainly not 
on the order of the equations) but on the regularity of the solution, including the 
regularity of the initial condition. 

2.3. Behavior on a parabolic case. We have performed a series of simple test 
on the linear parabolic equation, provided with periodic boundary conditions 



in the case where the viscosity /z is positive yet very small /i = 10~ 10 . We have chosen a 
standard spectral Fourier approximation in space (truncated series with N = 256 and 
chosen a second order in time implicit Euler scheme for a propagation over [0, T = 2]. 
The parareal in time algorithm is implemented with AT = 2.10~ 2 , dT = 10 -2 and 



The simulation are done in order to illustrate the behavior of the scheme as 
regards the regularity of the solution. We have thus chosen different types of initial 
conditions with different regularities by choosing uq(x) — Yl U(,e lix with 



In the case of low regularity: p — 0.5 or p = 1, the relative error after 15 parareal 
iterations is about 100! There is no blow up, the iterations can continue, the error 
eventually going to zero. For p = 4 or even better p — 10, the convergence is achieved 
after very few iterations, as is expected from Theorem 12.11 

These simple computations provide a preliminary illustration that the behavior 
of the parareal in time algorithm is deeply linked to the regularity of the solution. 

3. The periodic wave equation. 

3.1. The basic Fourier discretization. We consider the following one dimen- 
sional wave equation with periodic boundary conditions: Find u such that 




0, 

u (x), 

u(x + 27T, t) 



(2.11) 



St = 10" 4 . 




if I > 

if e = o 

if £ < 0, 




0, 

f(x), ^(x,0)=g(x), 

u(x + T,t), ^( X ,t) = ^(x + T,t) 



(3.1) 
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where / and g are periodic given initial condition, c stands for the speed of the wave, 
and T g M. + represents the period. 

The weak form of equation (|3.1[) is as follows: Find u such that \/t > 0,u(.,t) G 
S = {v G H\0,T),v(0) = v(T)} 

here, (u, u) = f Q uvdx. 

We use now the Fourier spectral method to semi discretize in space problem (|3.2p . 
That is, we choose 

S N = span{e u ^ x ,£=-N > --' ,N}, (3.3) 
and approximate u(x,t) by 

N 



UN 



(x,t)= J2 Mt)e U ^ x - (3.4) 



Note that ut(t) — ^(uN(x,t),e ie ^r~ x ), as results from the equality (e lf T x ; e lk T- x } = 
T(^,fe where J^fc denotes the standard Kroneker symbol. 

The solution un is defined by a Galerkin approximation so as to satisfy 

d 2 u N 2f du N dv N 

As is classical, this Galerkin approximation reduces to a system of differential equa- 
tions in the Fourier coefficients Uk(t). More precisely, let us choose iteratively vjy — 
e lfe T x ) k G [— AT, N] in (|3.5p . We then obtain the following system of independent 
equations 

T^(t) + ^Wu,(t) = 0, We [-JV, JV], (3.6) 

that can be exactely solved. After resolution, we know all the coefficients ui(t) in a 
closed form, then from equation (|3.4I) we get the Galerkin approximation of u(x,t) in 
Sn- 

Next, using respectively v = ^ in p.2[) and vjv = in (|3.5[) allows to derive, 
for any time t 

ll^(-,*)lli 2( o, T ) + c 2 H§^(-^)lli-(o, T ) = ll^(-,0)||i 2(0 , T) ^ c 2|||^ ( ., )||i 2(0 , ir) (3.7) 
and 

(3.8) 

respectively. If we now introduce the following Hamiltonian 

H(u N )(t)^\\^(.,t)\\l 2{0J) +c^^(.Ml H0 j) (3-9) 
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The conservation (|3.8p reads for any time t 

H{u N )(t) = H(u N )(0) (3.10) 

Of course the numerical discretization in time should preserve this invariant . 

Remark 2. As is often the case for Hamiltonian system, this basic invariant 
(which is here the total energy for the wave equation) is not the only one available, 
for instance it is immediately seen on \3. 6}) that by defining 

^W = l|l 2 + (f)WH 2 , (3.ii) 

for any £, £ = —N, .., N, we have 

Vt, h e (u e )(t) = h e (u e )(0) (3.12) 

which provides 2N + 1 invariant quantities (we notice that H = X)f=-jv ^i)- 

3.2. Illustration of the stability and instabilities of the plain parareal 
in time algorithm. In all what follows, the discrete solutions will always, for any 
time t, belong to Sn and will result from a Galerkin process. We will thus denote by 
u n the approximation (in time) of the value of uat(T„) (hence the full approximation 
of u(T n )). In order to test the parareal in time algorithm on this simple second 
order equation, we further discretize the semi discrete problem (|3.5j) by adding a time 
discretization, that is based, both for the coarse Gat and the fine J- at propagators, 
on the velocity Verlet scheme which is a second order in time, stable and consistent 
with the preservation of invariants (the wave equation being an Hamiltonian system, 
the symplectic nature of the velocity Verlet scheme makes it a method of choice) . 

The plain parareal scheme reads: 

41+1 = SaT«) foioN 
= 0At(< +1 )+.F AT (<)-<?AxK). ^ 

We perform simulations with the discretization parameters equal to 

1. T = 100, AT = 2.10- 1 , dT = 4.10" 3 , St = 2.10~ 6 . 

2. N = 30. 

for different initial conditions of different regularities: u (x) — ^ i £ Z U£e li ^~ x with 




if £ ^ 
if t = 0, 



with p p = 1, 5 or 10. The cases with p = 1 and p = 5 diverges while the case 
with p = 10 converges. 

3.3. Results with the plain parareal in time algorithm on a realistic 
test problem. In this example, the excitation is provided by a Ricker pulse, which 
is analytically defined as follows 

1. T = 5000 m, c = 2000 m/s 

2. and the initial conditions are f(x) = 0, g(x) = (l—2(f s n x ~ X3 ) 2 ) exp(— {f s -K^^ 
Here, f s = 2.5 Hz, is the frequency and x s — 2500 m is the position of the vibration 
source. 

The numerical results are obtained with a coarse propagator Q and fine propagator 
T that differ only with time step dT and 5t, respectively : 
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1. T = 100, AT = 2.HT 1 , dT = 4.10~ 3 , St = 2.10~ 6 . 

2. N = 30. 

The plot showing the evolution in time of the solution obtained after 15 iterations 
at different positions in the periodic space [0, ,T] on figure |3~T1 reveals an instability 
visible starting at time t = 80. The situation is even worse after 25 iterations. Actu- 
ally this is one case that illustrates the bad behavior of the parareal in time algorithm 
for this hyperbolic equation. Note however that for many other choices of the dis- 
cretization parameters 5t or dT or even N but the very same initial condition, the 
simulation is perfectly fine, stable and convergent. The instability is thus not system- 
atic, and changing slightly the discretization parameters, is often enough to make the 
scheme work. The plain parareal in time algorithm does lack of robustness, and this 
is what we want to correct in what follows. 




The exact solution shown on the next plot (figure [3~2l) illustrates the differences. 

The plots of the error between the exact solution and the solution obtained after k 
iterations, with k = 1, .., 15 (see figures l3~3l and l3~4")) illustrate clearly that the solution 
obtained by the parareal in time algorithm 

1. converges for short time, e.g. only 7 iterations are sufficient to recover the 
solution obtained with the sole sequential fine propagator at time equal to 
10. 

2. does not converge for longer propagations, here for a time equal to 100 

3. During the first 5 to 6 iterations the error of the parareal in time algorithm 
is dominated by the error of the coarse scheme, but after k = 6, the error 
increases largely to values that are order 1 and bigger 

4. After 15 iterations of the parareal in time algorithm the parareal solution is 
computed with an acceptable error until time t ~ 30 only. 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS 
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Fig. 3.2. Vibration at different positions obtained by the fine sequential method(8t = 2.10 - ®) 



wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS 



1 

0.1 
0.01 
0.001 
0.0001 
1e-05 
1e-06 
1e-07 
1e-08 
1e-09 
1e-10 




0.1 



10 



100 



Fig. 3.3. Error on the solution obtained by the parareal method(St = 2.10 e ,dT 
4.10 -3 , AT = 2.10" 1 ; 



Taking into account the variants of the parareal in time algorithm introduced in 
[7] to simulate Hamiltonian systems we infer that the bad behavior of the parareal 
scheme comes from the non preservation of the invariant quantities. Let us for example 
observe the evolution of the value of the quantity H defined in (|3.9I) 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS 
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Fig. 3.4. Error on the solution obtained by the parareal method(5t = 2.10 6 ,dT 
4.10 -3 , AT = 2.10" 1 ; 



The figures 13.51 and 13.61 illustrate the lack of conservation, that is coherent with 
the divergence of the scheme in this case. 

Finally we want to notice that we have performed the same simulation with a 
Ricker pulse of varying force (by modifying the value of / s ). The instability starts 
much earlier for larger f s confirming the relation between the instability and the 
singularity of the solution. 

3.4. Parareal algorithm with projection. Let us define the manifold M 
corresponding to the set of all functions v such that H(v) — H(u(0)). The Standard 
Projection Method on Ai (see chapter 4 of [13]) can be written as follows: For a given 
u, solve the system in (u, A) 

u = u + XS7H(u), 
where A € R is such that H(u) = H(u(0)), i.e. 

H(u + \VH(u)) = H(u(0)). (3.14) 

We then set 

u = ttm(u) (3.15) 

In the general case, equation (|3.14[) can only be solved by some iterative method, 
for example, Newton iterations, but in the current case, using the fact that H is 
quadratic, we notice and take into account that equation Q3.14[) is a second order 
polynomial in A that can be solved exactly. 

Then, combining the plain parareal scheme (|3.13[) with the projection method 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS 
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Fig. 3.5. Error on the energy obtained by the parareal method(St = 2.10 6 ,dT = 4.10 3 , AT : 
2.1Q- 1 ) 

wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS 




time 



Fig. 3.6. Error on the energy obtained by the parareal method(6t = 2.10 8 ,dT = 4.10 3 , AT : 
2.1Q- 1 ) 



above, we obtain the parareal method with projection to M: 



"n+1 
f.k+1 
"n+1 

"n+1 



= GatK) 

= ^at(<+ 1 ) + Jat«) - GatK), 



(3.16) 
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We now show the results obtained by parareal method with projection. We first 
note on figure [3~71 that the global solution after k = 15 iterations does not seem to 
present the instabilities of the plain parareal in time algorithm. The same conclusion 
holds after 25 iterations. 



wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1, PS + ProjH 
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Fig. 3.7. Vibration at different positions obtained by the parareal method with projection to the 
manifold M (St = 2.10" 6 , dT = 4.10 -3 , AT = 2.10 _1 J 



This is confirmed by the figures I3.8I and I3.9I where we note an improvement in the 
error between the exact solution and the discrete one as the iterations evolve. The 
error is always at worse of the same order as the error of the coarse propagator. 

Of course this is far from being what is expected, but the conservation of H 
does improve the behavior of the algorithm. The following figures I3.10I and I3.11I 
illustrate the good preservation of H along the trajectory for the parareal method 
with projection. Although the solution does not diverge as in the plain parareal case, 
it does not converge as we would like. 

In the results above, when doing the projection, we view all the Fourier modes 
as a whole system, and treat them altogether, with a unique Lagrange multiplier. 
Remember that each Fourier coefficient individually has an energy that is preserved 
as is explained in Remark [21 Instead of projecting on the manifold A4, where only one 
invariant is taken into account, we can have a manifold defined with more invariant 
properties. We have defined only two invariants based on two groups in the Fourier 
coefficients, those corresponding to high wave numbers and those corresponding to 
low wave numbers. That is, we partition all the Fourier modes into two groups, one 
contains the modes \k\ > K , the other one contains all the modes |fc| < K , we 
then perform a projection for each group separately to make the energy of each group 
conserved. 

The motivation comes from the fact that, when we discretize the solution by 
Fourier series, the contribution of the high frequency modes is much smaller than 
that of the low frequency modes. If we use the same Lagrange multiplier A to correct 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH 




time 



Fig. 3.8. Error on the solution obtained by the parareal method with projection to the manifold 
M (St = 2.10~ 6 , dT = 4.10 -3 , AT = 2.1CT ^ 



wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH 
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Fig. 3.9. Error on the solution obtained by parareal method with projection to the manifold M 
(St = 2.10" 6 , dT = 4.10 -3 , AT = 2.1CT 1 j 



all the modes, the contribution of the high frequency modes may be to magnify them 
and pollute the solution. 

So, let us introduce two Lagrange multipliers Ai and A2, one is for the high 
frequency modes, the other is for the low frequency modes only. 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH 
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Fig. 3.10. Error on the energy obtained by parareal method with projection to the manifold M 
(St = 2.1CT 6 ,dT = 4.1CT 3 , AT = 2.1(r 1 j 



wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH 
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Fig. 3.11. Error on energy obtained by parareal method with projection to the manifold M 
(St = 2.10" 6 , dT = 4.10 -3 , AT = 2.1CT 1 j 



The method is not more complex to implement. The numerical simulation still 
with the same N = 30 is based on K = 20; the results are much improved as can be 
seen on figures 13.121 13.131 and 13.141 

The solution itself, after 15 iterations, is in excellent accordance with the very 



16 



X. Dai and Y. Maday 





1e-10 




Fig. 3.13. Error on the solution obtained by the parareal with projection on two groups (St = 
2.10- 6 ,dT = 4.10 -3 , AT = 2.1Q- 1 ) 



fine approximation as can be seen on figure 13.151 

For the linear second order wave equation, the preservation of the invariant quan- 
tities based on the decomposition of the coefficients representing the solution into two 
groups: the low modes and high frequency modes, is thus a sufficiently mild ingredi- 
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wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH, 2 groups 
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Fig. 3.14. Error on the solution obtained by the parareal with projection on two groups (St 
2.1CT 6 , dT = 4.10 -3 , AT = 2.1Q- 1 ) 



wave equation: dt=2.0e-6, dT=4.0e-3, DT=2.0e-1 , PS + ProjH, 2 groups 
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Fig. 3.15. Vibration at different positions obtained by parareal method with projection on two 
groups (St = 2.10" 6 , dT = 4.10" 3 , AT = 2.10 -1 j 



ent to improve the behavior of the parareal in time algorithm. In the next section we 
shall extend this variant of the parareal in time algorithm to a nonlinear PDE. 

Note that, if when the method converges, i.e. Vn = 1, .., N, u\ — !• it^°, this limit 
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solution satisfies thanks to (|3.16[) 

that, here, due to the fact that the propagator J- at is symplectic and thus preserves 
the Hamiltonian, reads 

u n+l = ^AT(C) 

that proves that the method converges towards the fine solution, as expected. 
4. Burgers equation. 

4.1. The total Fourier/time splitting discretization. Let us consider in 
this section the following periodic viscous Burgers' equation in one dimension 

du d 2 u 1 d(u 2 ) „ . t 1N 

u{x, 0) = uq(x), 
u(x + T, t) = u{x, i), 

The viscosity v > is planned to be very small so that the solution will exhibit very 
sharp gradients. 

Taking into account that the periodic frame of the problem allows to use Fourier 
spectral approximation, we thus consider a semi discretization based on the Galcrkin 
Fourier approximation in Sn 

du N du N dv N 1 d(u N ) 2 

(-gptar) + K-g^-i ^) + ^-di-^ = °> WVN £ Sn - (42) 

We then choose to use a second order symmetrized splitting method based on the 
splitting between the diffusion effects (treated in the Fourier space exactly) and the 
nonlinear convection effects (treated in the physical space). The only difference be- 
tween the coarse propagator Q and fine propagator J- are relative to the time steps 
dT and fine time step 5t, respectively. Denoting by S either dT or St, the splitting 
schemes at each iteration proceeds as follows: 

1. Solve the following problem in Fourier space in the interval [t„ = n5, T n . i = 

*"»+§] 

du N du N dv N 

U N {r n )=u N ,n, (4-4) 

and denote u Nn+ i = un{t h+ i). Note that, as in the discretization of the 

wave equation, by choosing iteratively vn = e lk ^~ x , k G [-N, N], this results 
in a system of differential equations in the Fourier coefficients. 

2. Solve the following problem in [r„,r n+ i] 

(^,v N ) + l(^,v N )=0,Vv N eS N (4.5) 

u N {T n ) = u Nj7l+ i, (4.6) 

and denote uat(t„ + i) as u Nn+ i. Note that here, the solution procedure 
cannot be exact since the problem is non linear: the discretization in time 
is based on a 3 rd order Runge Kutta method and an exact evaluation of the 
Fourier coefficients of a< -^ JV - > by a Discrete Fourier Method based on S^jv- 
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3. Solve the following problem in Fourier space in [r n+ i , T n+ {\ 



du N du 



JV 



dv 



•n 



. )=0,Vv N eS N 
Ox Ox 

un{t ji+ i) = u Nn+ i. 



(4.7) 
(4.8) 



4.2. Results for the plain parareal in time algorithm. We have performed 
3 tests for the plain algorithm. First we note that, provided that the viscosity is large 
enough, the parareal in time algorithm is working fine. Indeed in case referred later 
as Case 1, we have chosen 

Case 1: 

1. fi = le - 2, 

2. iV c = 2 8 ,AT f = 2 8 



3. T = 2, AT = 2.10" 2 , cLT = 1.10" 



,8t 



1.10' 



The evolution is presented on figure 14.11 the relative error in the solution with 
respect to a solution computed with a much smaller time step is presented on figure 
4.21 showing that after 3 parareal iterations, the same convergence as for the fine 
sequential simulation is achieved. 



burgers: mu=1e-2, N f = 2 , T=2, DT=2e-2, dT=1e-2, dt=1e-4, fine 



= 




Fig. 4.1. Case 1: solution obtained by the parareal method(5t = le — A,dT = le — 2,DT = 2e — 2) 



For the second test, we have diminished the value of the viscosity. This case, 
referred to as Case 2, corresponds to the choices 
Case 2: 

1. fi = le — 3; 

2. N c = 2 9 ,N f = 2 9 , 

3. T = 2, AT = 2.10~ 2 , dT = 1.10" 2 , St = 1.10~ 5 
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1e-15 




Fig. 4.3. Case 2-2: relative error of the solution obtained by the parareal method(5t = le ■ 
5, dT = le - 2, AT = 2e - 2) 



For this case, the plain parareal method is not stable, as can be seen on the figure 
I4.31 after 3 iterations instabilities occur soon after time 1 where the solution of the 
inviscid Burgers equation becomes discontinuous. 
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It should be noted that, as explained in |17) . that using a coarser spacial repre- 
sentation helps in the stabilization of the process. Indeed, Case 3 below provides a 
stable parareal approximation 

Case 3: 

1. fi = le — 3; 

2. N c = 2 s , N f = 2 9 , 

3. T = 2, AT = 2.10" 2 , dT = 1.10" 2 , St = 1.1Q" 5 . 

As can be seen on the series of figures 14.41 14.51 and 14. 61 the method is stable, 
converges but very slowly. 

Burgers: mu=1e-3, N c = 2 8 , N f = 2 9 , dt=1.0e-5, dT=1.0e-2, DT=2.0e-2, PS 
1 , , , , , , , 




1e-12 - 



1e-14 1 1 1 — 1 — 1 1 — 1 — 1 1 — 1 — 1 

0.01 0.1 1 10 

time 

Fig. 4.4. Case 2-1: relative error of the solution obtained by the parareal method(8t = le — 
5, dT = le — 2, DT = 2e — 2) 

In the next sections, we shall improve the two last simulations by extending the 
method presented in the previous section. 

4.3. Results for the plain parareal in time algorithm with projection 
alone. To start with, let us notice that, due to dissipation, there is no quantity that 
is preserved along a trajectory, nevertheless, by using «jv as a test function in (|4.2[) , 
we derive that 

lkiv(t)||| 2( o >T )+2^^[^] 2 = ll"iv(0)||i 2(0;T) , 
and actually for any two times T n+ \ > T n > 

\\u N {T n+1 )f W) +2v ^ " +1 [^] 2 = IIMT n )||i W) . (4.9) 

In the parareal in time algorithm, starting from (remember that this represents 
the Galerkin approximation that should be denoted as v!^^ ) , we expect that 

IKUllMcT) - II^ATk +1 ]||L 2( 0,T) (4.10) 
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Burgers: mu=1e-3, N c = 2 B , N f = 2 a , dt=1.0e-5, dT=1.0e-2, DT=2.0e-2, PS 
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Fig. 4.5. Case 2-1: relative error of the solution obtained by the parareal method(5t = le 
5, dT = le — 2, DT = 2e — 2) 



Burgers: mu=1e-3, N c = 2°, N f = ¥, dt=1.0e-5, dT=1.0e-2, DT=2.0e-2, PS 
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Fig. 4.6. Case 2-1: relative error of the solution obtained by the parareal method(8t = le 
5, dT = le — 2, DT = 2e — 2) 



Since we do not want to compute Jat["' +1 ] at this stage, we approximate ||.Fat[w^ +1 ] |l 2 (o t) 

by 

1177 r„fc+llll ||-^AT[^n]||^(0,T) n„,fc + l|| , A „s 

KAT V U 2 (o,t) - — 17— Fn u « IU 2 (o,t)- (4-11) 

K U 2 (0,T) 



Stable parareal in time algorithm for hyperbolic system 



23 



What can thus be done is to consider the manifold corresponding to the 

set of all functions v such that |M| L 2(o,T) = ^pp^jf^ 2 ^ IU 2 (o,t) 

Second, we have noticed in the analysis of the wave equation that treating all 
the frequencies with a unique Lagrange multiplier is not effective enough. We thus 
partition the frequency space of functions in Sn into two (or more if necessary) groups 
Gi, i = 1,2 (or i = 1 , ..,1) defined by 

d = span{e u ^ x , N~ < \£\ < N+}, (4.12) 

where Lj" = 0, < Lf = L% - 1, = N + 1 (or Lj" = 0, < Lf = L% - 1, 
< L~\ = L^ +1 + 1, = N + 1 if more groups are around). The L 2 -projection on 
group Gi is denoted as Pi. In this case we consider the manifold [■M J ]„ +1 defined as 
being the intersection 

[^ +1 = iiH» +1 (4-13) 
where the group manifolds are defined as follows . 

r nfc+1 r o II r> I Ml [^t[\]] IU 2 (0,T) m , ,, , 

[m Tn = {v e S N , Pi («) U 2 (o,T) = rr^n — IK + U 2 (o,t)} (4-14) 

IKIU=(0,T) 

To any given v S Sjv, we associate the element Ti^nM-i^) defined by the Stan- 
dard Projection Method on [yVf 7 ]^" 1 " 1 , that involves the determination of / Lagrange 
multipliers. 

With this manifold the parareal scheme with projection reads 



Sat«) 

Sat« +1 ) + -Fat(u*) - Gat«), (4.15) 

/ ~ fc-|-l \ 

The results are as follows (see figures I4T71 14. 8[H3|) . We notice that the parareal in 
time algorithm is stabilized, but it does not converge much. Increasing the number of 
groups does not improve the behavior of this variant of the parareal in time algorithm. 

Remark 3. If we would know the exact value of the energy of each groups, the 
projection method would converge. The problem of course is that due to dissipation, 
and also to nonlinear effect that mixes up the low and high frequencies, the value of 
the exact energy in each groups is not known. This is the reason why we get into these 
difficulties. 

4.4. Results for the plain parareal in time algorithm with projection 
and damping. We notice that, for this non linear problem, the projection allows 
to stabilize the parareal in time algorithm but a little more is required. We have 
chosen to damp the high frequencies at each coarse step before proceeding to the next 
parareal iteration. This means that we introduce a damping function $ that, to any 
element vn in Sn provides an element wn — &{vn) in Sn, such that 

Wed, w t =vt (4.16) 
WeGi,i = 2,..I, We = f(i)v e (4.17) 

where f(i) is chosen such that high frequency parts are more and more filtered. In all 
the following results, we choose f(i) = jtti q > except for 2 groups, were we have chosen 



r <+i - 

\ "n+1 — 

u k+1 - 
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Burgers: mu=1 .Oe-3, N c =2 9 , N f =2 9 , dt=1 .Oe-5, dT=1 .0e-2, DT=2.0e-2, PS + ProjH, 2 groups 
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Fig. 4.7. Case 2-2: relative error of the solution obtained by the method parareal with projection 
to 2 groups(8t = le - 5, dT = le - 2, DT = 2e - 2) 



Burgers: mu=1 .Oe-3, N c =2 , N f =2 , dt=1 .Oe-5, dT=1 .0e-2, DT=2.0e-2, PS + ProjH, 2 groups 
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Fig. 4.8. Case 2-2: relative error of the solution obtained by the method parareal with projection 
to 2 groups(8t = le - 5, dT = le - 2, DT = 2e - 2) 



f(i) = 3 Q . With the same definition of the multigroup manifold the parareal 
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Burgers: mu=1 .Oe-3, N c =2 3 , N f =2 3 , dt=1 .Oe-5, dT=1 .0e-2, DT=2.0e-2, PS + ProjH, 2 groups 
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Fig. 4.9. Case 2-2: relative error of the solution obtained by the method parareal with projection 
to 2 groups(8t = le - 5, dT = le - 2, DT = 2e - 2) 



scheme with projection and damping reads 



"n+1 



u 



5atK) 

<?AT« +1 ) + FMO - QAT(u k n ), 



V k+1 - IT L+l(il k + 1 ) ( 4 ' 18 ) 

"n+1 — [M'\n ^ n+l) 

<X\ = 

The following figures I4.10[ 14.111 and 14.111 present the results of this improved 
parareal method. We notice that the parareal method is now converging. 

Remark 4. Note that the strategy used here consisting in adding a little bit of 
dissipation on the high modes to the parareal scheme goes in the direction suggested 
in l^/ and is coherent with the spectral viscosity proposed in \22^ . Note also that the 
solution should be measured after the projection and before the damping. Note finally 
that the fine propagator acts on the undamped solution so that the dissipativity that 
is added does not affect the quality of the solution (e.g. as any parareal scheme, after 
N th iterations, the solution is the same as the one of the plain fine propagator) but 
improves the quality of the scheme. 

The projection with more than two groups performs slightly better, see e.g. figure 
14.131 where only the results for iterations between 7 and 15 are presented. It may be 
interesting to keep this in mind even though the separation in more than two groups is 
not completely straightforward in the non periodic case where finite difference, finite 
differences or finite volume techniques will be used. 

5. Preliminary analysis. In this section, we provide some simple elements 
of analysis that prove the stabilty of the previous approximation with the parareal 
method with projection (and damping). 

Let us first assume that the fine propagator is L 2 stable in the sense that there 
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Burgers: mu=1 .Oe-3, N c =2 , N f =2 , dt=1 .Oe-5, dT=1 .Oe-2, DT=2.0e-2, PS + ProjH + Damp, 2 groups 
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Fig. 4.10. Case 2-2: relative error of the solution obtained by the method parareal with projec- 
tion to 2 groups with damping(5t = le — 5, dT = le — 2, DT = 2e — 2) 



Burgers: mu=1 .Oe-3, N c =2 , N,=2 , dt=1 .Oe-5, dT=1 .Oe-2, DT=2.0e-2, PS + ProjH + Damp, 2 groups 
1 




Fig. 4.11. Case 2-2: relative error of the solution obtained by the method parareal with projec- 
tion to 2 groups with damping(8t = le — 5, dT = le — 2, DT = 2e — 2) 



exists a constant a > independent of the discretization parameters such that, 

WeZ 2 (0,T), ||J r AT(«)|U» ( o,T)<(l + aAr)||»|| £3{0iT) (5.1) 
Note that both for the wave equation and the Burgers' equation, the previous 
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Burgers: mu=1 .Oe-3, N c =2 , N f =2 , dt=1 .Oe-5, dT=1 .Oe-2, DT=2.0e-2, PS + ProjH + Damp, 2 groups 



* 0.0001 




1e-12 



1e-14 



Fig. 4.12. Case 2-2: relative error of the solution obtained by the method parareal with projec- 
tion to 2 groups with damping(5t = le — 5, dT = le — 2, DT = 2e — 2) 



Burgers: mu=1 .Oe-3, N c =2 , N,=2 , dt=1 .Oe-5, dT=1 .0e-2, DT=2.0e-2, PS + ProjH + Damp, 1 groups 
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Fig. 4.13. Case 2-2: relative error of the solution obtained by the method parareal with projec- 
tion to 10 groups with damping(St = le — 5, dT = le — 2, DT = 2e — 2) 



inequality holds with a — 0. 

Let us start with the case where no dumping is implemented, then, from (|4.15|) 
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Burgers: mu=1 .Oe-3, N c =2 a , N f =2 a , dt=1 .Oe-5, dT=1 .Oe-2, DT=2.0e-2, PS + ProjH + Damp, 1 groups 
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Fig. 4.14. Case 2-2: relative error of the solution obtained by the method parareal with projec- 
tion to 10 groups with damping(5t = le — 5, dT = le — 2, DT = 2e — 2) 



we derive 



\u k+1 



+1 — 7T [M I ]n +1 Hi 2 (0,T) — 2J ll-^[ M n+l]lll 2 (0,T) 

2 

L 2 (0,T) ||„,fc + l||2 



^I!^at[<]]" 2 



E 



fcTTl II % llL 2 (0,T) 



i=1 ll u 'nllL2( 0jT ) 

_ ij^fep H u ™ Hl 2 (o,t) 

H"nllL 2 (0.T) 

<(l + aAr) 2 ||^ +1 || 2 L2(0 , T) (5.2) 

This proves the stability by induction 

ll<illU 2 (o,T) < (1 + aAT) n+1 ||uq +1 = uo||l 2 (o,t) 

< e QT KIU 2 (o,T) (5.3) 

6. Concluding remarks. In this paper we have analyzed and proposed a cure to 
the lack of robustness of the parareal in time algorithm applied to hyperbolic systems 
or convection diffusion problems with small diffusion. The lack of robustness, leading 
sometimes to instabilities of the algorithm before convergence is a consequence of the 
lack of regularity that the solution develops when time runs. The cure is inherited 
from the extension of the parareal in time algorithm for Hamiltonian system and 
consists in projecting the solution over an energy manifold defined on the fly, and 
improved iteratively, where the exact solution lies. This simple procedure allows to 
stabilizes the algorithm that now converges after some iterations. It should be noted 
that, similarly as for other types of equations, the convergence of the parareal in time 
algorithm allows to speed up the solution procedure but is not optimal in terms of 
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parallel efficiency since, roughly speaking, a convergence that is achieved after 10 or 15 
iterations leads to a parallel efficientcy of A or at best, which is far from optimal. 
It is not the purpose of this paper to provide a full optimal implementation of the 
approach. Let us just notice that the combination of the parareal in time algorithm 
with standard domain decomposition methods is a way to improve the efficiency and 
use more of the processors that are available on current architectures than the one 
that a plain domain decomposition technique is able to optimally use. We refer to 
[12"] for more on this matter. 

Acknowledgments. This work has been supported in part by the Agence Na- 
tional de la Recherche, under the grant ANR-06-CIS6-007 (PITAC). 



REFERENCES 

[1] L. Baffico, S. Bernard, Y. Maday, G. Turinici, G. Zerah, Parallel in time molecular 
dynamics simulations, Physical Review. E 66(2002) 

[2] G. Bal , On the convergence and the stability of the parareal algorithm to solve partial differen- 
tial eguations, In Domain Decomposition Methods in Science and Engineering, Kornhuber 
R, Hoppe RHW, Keyes DE, Periaux J, Pironneau O, Xu J (eds). Lecture Notes in Com- 
puter Science and Engineering, vol. 40. Springer: Berlin, 2004; 425432. 

[3] G. Bal, Y. Maday, A "parareal" time discretization for non-linear PDE's with application 
to the pricing of an American put, in Recent Developments in domain Decomposition 
Methods, L. F. Pavarino, A. Toselli (eds), Lecture Notes in Computational Science and 
Engineering, 23, Springer, 189—202 (2002) 

[4] K. BuRRAGE, Parallel and seguential methods for ordinary differential equations, Numerical 
Mathematics and Scientific Computation, Oxford Science Publications, The Clarendon 
Press, Oxford University Press, New York, 1995. 

[5] K. BuRRAGE, Parallel methods for ODEs, Advances in Computational Mathematics 7(1-2), 
1-3, 1997. 

[6] F. Chouly and M. A. Fernandez, An Enhanced parareal Algorithm for Partitioned Parabolic- 
Hyperbolic Coupling, AIP Conference Proceedings, 1168 1, pp 1517-1520, 2009. 

[7] X. Dai, C. Le Bris, F. Legoll, Y. Maday , parareal algorithms for Hamiltonian systems, 
submitted to M 2 AN 

[8] C. Farhat, M. ChandesriS , Time- decomposed parallel time-integratorspart I: theory and 
feasibility studies for fluid, structure, and fluidstructure applications, International Journal 
for Numerical Methods in Engineering 2003; 58(9): 13971434. 
[9] C. Farhat, J. Cortial, C. Dastillung, H. Bavestrello , Time-parallel implicit integra- 
tors for the near-real-time prediction of linear structural dynamic responses, International 
Journal for Numerical Methods in Engineering 2006; 67(5):697724. 

[10] J. CORTIAL, C. Farhat , A time-parallel implicit method for accelerating the solution of non- 
linear structural dynamics problemslnternational, Journal for Numerical Methods in En- 
gineering 2009; 77 -.451470. 

[11] P. F. Fischer, F. Hecht, and Y. Maday. , A parareal in time semi-implicit approximation of 
the navier-stokes equations, In Domain decomposition methods in science and engineering, 
volume 40 of Lect. Notes in Comput. Sci. Eng., pages 433440. Springer, Berlin, 2005. 

[12] R. Guetat. , Mthode de parallelisation en temps : Application aux methodes de decomposition 
de domaine PhD Thesis, UPMC Paris 2011. 

[13] E. Hairer, C. Lubich, and G. Wanner, Geometric numerical integration: structure- 
preserving algorithms for ordinary differential equations, ?? 

[14] J.-L. Lions, Y. Maday, G. Turinici , Resolution d'EDP par un schema en temps "parareel" , 
Comptes Rendus de l'Academie des Sciences de Paris — Serie I — Mathematique 2001; 
332(7) :661— 668. 

[15] Y. Liu, and J. Hu„ Modified propagators of parareal in time algorithm and application 
to Princeton Ocean model, International Journal for Numerical Methods in Fluids, 57: 
17931804. (2008). 

[16] Y. Maday. , parareal in time algorithm for kinetic systems based on model reduction, In High 
dimensional partial differential equations in science and engineering, 41 CRM Proc. Lecture 
Notes, pages 183 — 194. Amer. Math. Soc, Providence, RI, 2007. 



30 



X. Dai and Y. Maday 



[17] Y. Maday. , The "parareal in time" algorithm, In Sub-Structuring Techniques and Domain 

Decomposition Methods, 24 Computational Science, Engineering Sz Technology Series, 

pages 19 — 44. Saxe-Coburg Publications, 2010. 
[18] Y. Maday, G. Turinici , A parareal in time procedure for the control of partial differential 

equations C.R.Acad Sci. Paris Ser. I Math 335 387—392 (2002) 
[19] D. Mercerat and L.Guillot and J. -P. Vilotte , Application of the parareal Algorithm for 

Acoustic Wave Propagation, AIP Conference Proceedings, 1168 1, pp 1521-1524, 2009. 
[20] J. Nievergelt. , Parallel methods for integrating ordinary differential equations, Commun. 

ACM, 7(12):731—733, 1964. 
[21] G. Staff and E. R0NQUIST , Stability of the parareal Algorithm In Domain decomposition 

methods in science and engineering, volume 40 of Lect. Notes in Comput. Sci. Eng., 

pages 449 — 456. Springer, Berlin, 2005. 
[22] E. Tadmor, Convergence of Spectral Methods for Nonlinear Conservation Laws, SIAM Journal 

on Numerical Analysis Vol. 26, No. 1 (Feb., 1989), pp. 30-44 



